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,__i Abstract 

o: 

C^l Treating a realistic problem in any field of reaction kinetics raises a series of problems: we review these illustrated with examples 
^ using ReactionKinetics, a Mathematica based package. 
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1. Introduction 



In Part I of our paper 



Toth et al 



|2011] we have formulated 



a: 



the requirements for a reaction kinetics package to be useful for 



a wide circle of users. 

In the present Part II we enumerate the major problems aris- 
ing when writing and using such a package. It turns out that the 



GO solution of some problems is far from being trivial. 

Ov 

O; 

. 2. Parsing 

o ■ 

o; 

As we mentioned in the first part of our paper, at the very 

. . beginning of application of computers to chemical kinetics the 

> ■ 

• ^h problem arouse how to construct the induced kinetic differen- 

%4 tial equation of a reaction without making to much errors. The 
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best thing is if the chemist provides the reaction steps, and the 
program creates the induced kinetic differential equations auto- 
matically, if the kinetics is supposed to be of the mass action 
type. If it is not then reaction rates should also be provided by 
the user. 

As our program is based on a modern mathematical program 
package (often called by the misnomer computer algebra sys- 
tem) there is no limitation as to the size of the reaction to be 
handled, see the Lendvay example in our previous paper. Con- 
fer this with the quite typical restriction ". . . has the capabilities 
for handling up to 590 differential equations." 

It may be instructive to cite the code building the right hand 
side of the induced kinetic differential equation of the reaction 
(note that reversibility is not assumed) 



M 



M 



Y a(m, r)X(m) -A ^j3(m, r)X(m) (r = 1, 2, . . . ,R) 



m=\ 



m=\ 



given the matrices (a,j8) of molecularities, the vector of reac- 
tion rate coefficients (k = (& r )* =1 ) and the concentration vector 
(c) the right-hand side of the kinetic ODE gets an exceedingly 
transparent form (in coding as well): 



08-a).(kTimes( 



c a ), 



where notice that both dot and componentwise products are 
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used. The reaction step vector y(-, r), where y — fi — a ex- 
presses as the effect of the r th reaction step. 

3. Combinatorics 

In reaction kinetics graphs of many types are useful and used, 
though it is not quite obvious how to visualize them to be 
the most informative for users, especially concerning to large 
mechanisms. Hence we accept the representations offered by 
Mathematica. 

The Volpert graph is one of the widely applied graphs in re- 
action kinetics which proved to be a good tool to investigate 
problems emerging e.g. in the theory of kinetic ordinary dif- 
ferential equations. It is a directed bipartite graph with species 
and reaction steps as its vertices and with a(m, r) edges from 
the vertex representing species X(m) into the vertex represent- 
ing reaction step r and with /3(m, r) edges from the vertex rep- 
resent ing reaction step r, i nto the vertex representing species 



X(m) [Erdi and Toth, 



1989]. At the same time this is the graph 



one can see in textbooks on biochemistry, in metabolic maps 
or at MaCKiE Workshops. With out any additional advan tage 



they can also be called Petri nets ITMartinez and S ilva. 



198211 . A 



model of glycolysis is a built-in reaction of our program pack- 
age therefore having got it we can easily draw its Volpert graph 
which turns out to be not so much different from the one used 
in textbooks on biochemistry (see Figure[3). 

4. Linear (Diophantine) equations 

Complex reactions occur via pathways of simple "elemen- 
tary" reaction steps. It is a common problem that the overall (or 
global) reaction is measured, and one would like to reconstruct 
the underlying network of simple reaction steps. By simple here 
one may mean reaction steps of order not h igher than two. Let 



us consider an example IKovacs et al 



2004]. We start form the 



overall description of the permanganate/oxalate reaction 

2Mn0 4 ~ + 6H + + 5C 2 H 2 4 ^ 2Mn 2+ + 8H 2 + 10CO 2 . 

The following steps lead to its the decompositions to elemen- 
tary steps: 



Volpert graph or Petri net 



1-I.dBP^ DHAP+GADP 



2PG^ H20 + PPP 



I3BPG 
1.3BPG-.GADP hi! .! 





/ 

2PG-.3PG j L 

1 

3PG 

\ 


3PG -> 2PG 




Figure 1 : The Volpert graph or Petri net of a glycolysis model 



H 2 C 2 4 HC 2 4 W ClO^ 



Mn 



TT 



MnC 2 4 



Mn0 4 Mn0 2 Mn 3+ C0 2 H 2 COj 
[Mn0 2 ,H 2 C 2 4 ] [Mn(C 2 4 )] + [Mn(C 2 4 ) 2 r 

[MnC 2 4 ,Mn0 4 ,H + ] [MnC 2 Of , MnO;] + 
[MnC 2 Of , MnO v H + ] 2+ [H + , MnQ 2 , H 2 C 2 Q 4 ] + 

Table 1: Species of the permanganate/oxalic acid reaction according to 



Kovacs et al 



200411 . Species with initially positive concentration are typeset 



bold. 



1 . determine the combinatorially possible species and select 
from those the chemically acceptable ones, 

2. determine the combinatorially possible reaction steps 
which are simple in the above sense and select from those 
the chemically acceptable ones, 

3. find those representations of the given overall reaction 
which are chemically acceptable. 

As we have seen in Part I of our paper, the second and third 
steps reduce to the solution of line ar Diophantine equations. 



The authors of IKovacset al. , 



2004] worked with 19 species, 



including four hypothetical transient species, made up of four 
atomic constituents (and charge); these are listed in Table Q] 
The corresponding atomic matrix (see part I) is a 5 x 19 matrix. 
To generate all combinatorially feasible elementary reactions, 
2 • 19 + (2) = 209 systems need to be solved. The results 
give 1022 combinatorially feasible elementary steps. Based 
on chemical evidence a large number of them were eliminated, 
leaving 673 elementary reactions. 

The same computat ional results were reproduced in 



I Pappand Vizvari . 



2006], where further combinatorial analy- 



sis revealed that only 297 of the previously generated elemen- 
tary reactions may be part of decompositions, and that 3 of the 
species cannot be generated during the reaction. This analysis 
takes into account which five species are present initially in the 
reaction (based on experiments), and is an adaptation of the in- 



dexing process by 



Volpert 1 1972]. It is perhaps worth mention- 



ing that the result (in this specific example) is not sensitive to 
the assumption on the initial species; the same elementary steps 
and species remain if every non-complex species is assumed to 



be present initially in the vessel. 

Decompositions of the overall reaction can be obtained via 
the solution of another linear Diophantine system of equations. 
With the above analysis the dimensions of this system are re- 
duced to 16 x 297. The generation of decompositions can be 
further accelerated by preprocessing: simple analysis, using lin- 
ear programming, shows that every decomposition consists of 
at least 15 steps, and that two elementary steps take part in ev- 
ery decomposition, one of them with a coefficient of at least 
four. After this preprocessing, it is possible to find all decom- 
positions of 15-17 steps, and thousands of more complex de- 
compositions, from which chemically acceptable ones can be 
selected. 

We refer the reader to the two papers cited above for details 
of these computations, and a review of a number of algorithms 
for the solution of linear Diophantine equations. All of these 
algorithms and the various methods of combinatorial analysis 
mentioned above are now built in to our program. 

5. Numerics 

5.1. Stiffness 

One of the major problems when solving initial value prob- 
lems describing chemical reactions is — as in general — stiffness: 
the situation when the reaction proceeds along more than one 
time scales, reaction rates of different steps are of different mag- 
nitude. Then the codes meet the dilemma: if they choose a 
small step size, then the numerical solution will need much 
time, if they take a too long step size, then interesting phe- 
nomena occurring in the time evo lution of th e fast species 



will be lost. Starting w ith Gear [Gear, 



I But cher and Cash , 



1992] and Butcher 



199011 many methods fulfilling the require- 



ments of the chemist have been elaborated and built in into pro- 
gram packages. The following example show that such an in- 
sertion has been carried out especially successfully in the case 
of Mathematica. 



The Robertson problem IIRobertson , 



1966] 



. 004 „ „„ 3xl0 7 „ „ „ „ 10 4 , „ 

A — > B, 2B — >B + C, B + C — > A + C 



Mass conserved 




Figure 2: The Volpert graph of the Robertson model 

is a popular benchmark problem of stiff type, because there is 
a nine order of magnitude difference between the reaction rate 
coefficients. To make matters worse, the solution is required on 
the time interval [0, 10 11 ]. Figure I5TTI shows its Volpert graph. 



total cone. 
1.000010 r 



200000 400000 600000 800000 1 x 10* 

Figure 3: The total mass in the Robertson model 

Finally, we remark that a symbolic approach to unveil stiff- 
ness when it does not manifest itself on the surface has been 
given by Pro f. Goldshtein at the conference, and see also 



I Bvkov et al. 



2008] 



No matter what the values of the reaction rate coefficients are, 
the derivatives of the concentrations sum up to zero, meaning 
that total mass is conserved. One way to verify this fact is to 
execute the following: 

Total [RightHandSide ["Robertson", ki , k 2 , k 3 ]], 

which gives 0, cf. [?]. The concentrations are calculated 
numerically as simply as usual (with our built-in function 
Concentrations): 

Concentrations [{"Robertson"}, J0.04, 3xl0 7 , 10 4 ), 
{1, 0, 0}, {0.10 11 }] 

Also the numerical method is good enough to conserve mass 
(Figure l5TTT l. 



a decreasing, c increasing 




^- log(time) 



Figure 4: Species A and C in the Robertson model 



Let us see the individual concentrations of species A and C us- 
ing the logarithmic time scale (Figure \5A\ . But what about the 
species B? See Figure I5TI 

One of our colleagues, Dr. R. Horvath, was so kind as to 
solve the same problem using M ATLAB; the same results were 
obtained within the same CPU time, but much less automati- 
cally. 



5.2. Bifurcations 

A usual method to show that oscillatory solutions exist in the 
induced kinetic differential equation of a reaction is to apply the 
well-known Andronov-Hopf bifurcation theorem. Visualizing 



The quantity of the intermedier is so small 
that we were unable to see it in the previous figures. 



7. Stochastic simulation 
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Figure 5 : Species B in the Robertson model 

the effect of the change of the parameters is really simple using 
the Manipulate command of Mathematica see also the last 
example in the first part of our paper. If you are interested how 
easily this can be done you might have a look at the demon- 



stration by IVardai and Tothl In fact the method used there is 
based on particular ideas and has not as yet been built in into 
our program package handling general models. 



6. Symbolic methods 



Symbolic solutions are easier to find nowadays, still it is 
a hard problem even in the case of stationary points, see 
our previous paper. The most interesting developments are 
expected in these fields, e.g. how to find automatically 
the equivalence (of some kind) of different kinetic models 



I Mendez and Femal . 



201111 . how to show the possibility or im- 



possibility of transforming o ne dynamics into anoth er by (say, 



orthogonal) transformations |T6th and Hais, 



1986a] using per- 



haps t he theory of algebraic i nvariants of differential equa 



tions IHalmschlager et al 



200411 . how to find a minimal model 



with given properties within a class of reactions [Wilhelm, 



2009; 



Toth and Hars, 



1986b], given an induced kinetic differ- 



ential equation how to find a reaction with a give n structure 



I Szederkenvi , 



2010; 



Szederkenvi and Hangos , 



In reaction kinetics, it has always been obvious from the the- 
oreti cal point of view that in the case of small systems (see 



e-g- 



Aranvi and Tothl B1977I0 . and systems operating around 



an unstable stationar y state (e.g. models of chirality, see e.g. 



Barabas et al. 



|2010]) are better described with the standard 



continuous time discrete state stochastic model of chemical re- 



actions, see e. g. [Erdi and Toth, 



19891 Chaptr 5]. However, 



it turned out to be much more relevant only recently as with 
the development of analytical techniques it became possible to 
measure individual molecules. That is the reason why simula- 
tion of the stochastic model forms a relevant part of our pro- 
gram package. 

Here we are considering only the usual continuous time dis- 
crete state model which is a Markovian jump process. Denote 
by X(f) the number of species present in the system at time t. 
Now, our investigated process is defined to evolve in time so 
that the law 

P (/?, takes place in [t, t + h) | X(0) = K r (X(f)) h + e(h)h ( 1 ) 

is fulfilled where e{h) — > if h — » and k/s (r - 1,2,...,/?) 
are a kind of combinatorial functions of the vector of num- 
bers of species similar to but slightly different from the prod- 
uct of powers (Kurtz kinetics) involving also the reaction rate 
coefficients. From this assumption one can easily determine 
the infinitesimal generator of this Markov process from which 
the Kolmogorov forward and backward equations as well as 
the master equation are obtained. However, our formulation 
is equivalent to the following explicit representation: 

x(t) = x(0) + £ z ' (f Kr(X(s)) ds ) r0 ' r) (2) 

where the Z/s are independent, unit-rate Poisson processes. 
Note that there are many known stochastic simulation algo- 
rithms for the presented problem, but we intend to summarize 



2011], etc. 



only a few of them, namely: the direct method dSipos et al 
1 1974allbl0 : let k(x) := £f = i K r (x) and suppose T > is given, 
then the initialization step is: x(0) := Xo e K M , t := then the 
algorithm works as follows: 



Whileff < T, 

y :=RandomChoice[{^^;r = 1,2, ...,/?} -> 

{R r ;r= 1,2,... ,/?}]; 

.? : = RandomVariate[ExponentialDistribution[^(x(f))]]; 

x(f + s) := x(f) + y; ? = t + s] 



Notice that the idea of the direct method relies on the def- 
inition ([TJ. Also there exists several variants (e.g. first reac- 
tion method) and improv ements of this early method (see e.g. 



Gibson and Br uck [2000]). In what follows we briefly present 



the approximation methods, which have their roots in formula 
©. At the heart of these kinds of methods is the leap con- 
dition, that is choose r to be small enough so that the change 
in the state in [t, t + t) causes no relevant change for the K r 's. 
Again, the initialization step is x(0) := Xq e M. M , t := then 
the algorithm works as follows: 
Whileff < T, 

Choose t so as to satisfy the leap condition; 
Do[p r := RandomVariate[PoissonDistribution[A: r (x(f))T], 
{r,l,J?}]; 
x(t + s) := x(t) + J* , J P,-(x(f), x(f + r), Pr y, 

t -t + T] 



• If we take P r (x(t),x(t + f),p r ) '■- p r , then we get the ex- 
plicit T-leaping method. 

• If P r (x(t),x(t + f),p r ) := Pr - TK r {x(t)) + K r (x(t + t)), We 

are led to the implicit T-leaping method. 

• If we set P r (x(t), x(t + t), p r ) := p r - |/c r (x(?)) + jK r (x(t + 
t)) then what we get is t he trapezoidal T - leapin g method 



dCao and Petzo ld [2005]; 



Rathinam et al. 



12003]). 



However, these methods may also be considered as the stochas- 
tic analogues of certain numerical schemes applied for ordinary 
differential equations. Choosing the appropriate t has turned 
out to be a non-trivial task. It is also a challenging problem 
to a void negative pop ulation for x(f) in the simulation process, 



to handle stiff problems, where typically implicit approaches 
give the most appropriate results. These and several other algo- 
rithms have also been implemented into our program package. 
We also included conversion of units functions which provide 
an elegant way to compare the solution of the induced kinetic 
differential equation with the results of the stochastic simula- 
tion process which is considered to be taking place in a certain 
volume. In Figure|7]we can see this comparison on the example 
of Brusselator model: 



A -^ X,X — » Y,2X + Y -^ 3X,X -^ P 



where A and P are external species. 



Species: X. Y 




Figure 6: Stochastic Brusselator model is considered with reaction rate coeffi- 
cients: ki = 1.92, £2 = 5.76, iq = 5.6, k$ = 4.8 and initial conditions .To = 500, 
yo = 720. The "noisy" curve describes the evolution of the stochastic model, 
where the volume is 10~ 21 dnr . 



In Figure [7] one can see the behaviour of the Autocatalator 
model: 



k\ kn ki k* 

A -U X,X -A Y,X + 2Y — ^ 3Y, Y -A P, 



where again A and P are external species. 



see flCao et al 



200511 . Especially these last methods are used 



number of species 



Species: X, Y 




Figure 7: Stochastic Autocatalator model with reaction rate coefficients: k\ = 
110.2, k 2 = 0.094, k-i = 0.01 1, U = 90.34, where the initial conditions jr = 15, 
j>0 = 80 were considered. 

8. Estimation methods 

8.1. Based on the deterministic model 

Estimating the reaction rate coefficients (or even the Arrhe- 
nius parameters fcn, n and A in the formula 

k(T) = k T n exp(-A/(/?T)), 



see IINagv and T. . 



2011]) based on measurements is a kind of 
art. Here we only make a few remarks. 

First of all, Mathematica gives the possibility to obtain esti- 
mates (seemingly) without any kind of iteration: it is capable 
using the numerical solution of a differential equation in the 
same way as an explicitly given function. Let us start from the 
deterministic model of the reaction 



Then, the model needed by the function FindFit 

FindFit [data, model [a, b] [x] , 

{{a, 0.7}, jb, 0.2}}, x] 

is defined in the following way. 

model [a_,b_] := (model [a, b] = First [x /. 

NDSolve[x' [t] == b x[t] - a x[t] 2 , x[0] == 2}, 

x, {t, end}]]) 

As a result we get (0.344129,0.75174) instead of 
(0.33, 0.72), not a bad result, but certainly this is only a toy 
example. The agreement between the "measurements" and the 
fitted model is quite good, see Figure [8TI A systematic inves- 
tigation of similar (and further) estimation procedures for more 
complicated models is in progress. 

No worse fitting ever! 




Figure 8: Fitted model and data 



2X ^ X, 



where k\ - 0.33 and k-i - 0.72 with the following initial 

concentration of X: x(Q) - 2. Let us generate experimental data 

from the numerical solution and add a small amount of error. 

end = 7; sol = First [ x /. 

NDSolve[{x' [t] == -0.33 x[t] 2 + 0.72 x[t], 

x[0] == 2}, x, {t, end}]] ; 

times = N [Range [0, 5 end]/ 5]; 

data = Transpose [{times, sol [times] + 

RandomReal[.01, 5 end + 1]}]; 



We also mention that a method similar to simulated anneal 
ing applied to the determination of Arrhenius para meters of re 



action rate coefficients has also been elaborated INagv et al 



20111. 



8.2. Based on the stochastic model 

Based on the stochastic model one can obtain estimates in 
different ways. If we are able to collect data on the individ- 
ual changes of the molecules what is not so impossible in these 
days, then we can have maximum likelihood estimates. What 
is more even equilibriu m fluctuations migh t be enough to esti- 



mate the parameters, lErdi and Toth , 



1989, Subsection 5.8.7]. 



Another approach is to u se a kind of implicit linear regression: 



Hangos and Toth 



1 1988]. Methods elaborated to estimate pa- 
rameters of mass communic ation can also be adapted to the 
goals of reaction kinetics, see[ 



Kovacs. 



9. Visualization 

No external package is needed to provide camera ready fig- 
ures for papers in any of the following forms: eps, jpg, pdf, 
bmp, gif, wmf etc. (see a delicate example on Figure|9j. 




Figure 9: The deterministic Brusselator model plotting in 2D using 
Manipulate in which one can change the initial conditions in real time (us- 
ing the "locator" object) 



10. Discussion, further plans 

We allow non mass action type kinetics even now, therefore 
it does not seem to be a hard task to include temperature de- 
pendence, so important in the case of modelling e.g. atmo- 
spheric chemistry and combustion. Surely, reaction diffusion 
equations should also be treated in full generality including 
symbolic treatment, as well. The reduction of the number of 



variab les is also a very important topics, see e.g. JToth et al 



Finally, please visit the website: 

demonstrations .wolfram, com and look for the demonstra- 
tions of Attila Laszlo Nagy and Judit Vardai. 

10.1. Acknowledgments 

The authors thank the cooperation with dr. R. Horvath. 
Many of the participants of MaCKiE 201 1 contributed with in- 
centive ideas. Further requirements, criticism and problems to 
be solved are wanted. 

References 

P. Aranyi and J. Toth. A full stochastic description of the Michaelis-Menten 
reaction for small systems. Acta Biochimica et Biophysica Hungarica, 12 
(4):375-388, 1977. 

B. Barabas, J. Toth, and G. Palyi. Stochastic aspects of asymmetric autocataly- 
sis and absolute asymmetric synthesis. The Journal of Mathematical Chem- 
istry, 48(2):457^189, 2010. 

J. C. Butcher and J. R. Cash. Towards efficient Runge-Kutta methods for stiff 
systems. SIAM Journal on Numerical 'Analysis, 27(3):753-761, 1990. 

V. I. Bykov, V. Gol'dshtein, and U. Maas. Simple global reduction technique 
based on decomposition approach. Combustion Theory and Modelling, 12 
(2):389-405, 2008. 

Y. Cao and L. Petzold. Trapezoidal r-leaping formula for the stochastic simula- 
tion of biochemical systems. Proceedings of Foundations of Systems Biology 
in Engineering, pages 149-152, 2005. 

Y. Cao, D. T Gillespie, and L. Petzold. Avoiding negative populations in ex- 
plicit poisson r-leaping. The Journal of Chemical Physics, 123(5), 2005. 

P. Erdi and J. Toth. Mathematical models of chemical reactions. Theory and 
applications of deterministic and stochastic models. Princeton University 
Press, 1989. 

C. W. Gear. Invariants and numerical methods for odes. Physica D, (60):303- 
310, 1992. 

M. A. Gibson and J. Bruck. Efficient exact stochastic simulation of chemical 
systems with many species and many channels. The Journal of Physical 
Chemistry A, 104(9):18761889, 2000. 

A. Halmschlager, L. Szenthe, and J. Toth. Invariants of kinetic dif- 
ferential equations. volume 1, pages 1-14, Szeged, 2004. URL 



19971 



http: //www. emis.de/ journals/EJQTDE/7/714. html Proc. 7'th 

Coll. Qualitative Theory of Diff. Equ. 
K. M. Hangos and J. Toth. Maximum likelihood estimation of reaction-rate 

constants. Computers and. Chemical Engineering, 12(2/3):135-139, 1988. 
B. Kovacs. An intensity estimation method for inhomogeneous point processes 

and its application in reaction kinetics and telecommunications. Submitted. 



K. Kovacs, B. Vizvari, M. Riedel, and J. Toth. Computer assisted study of the 
mechanism of the permanganate/oxalic acid reaction. Physical Chemistry 
Chemical Physics, 6(6): 1236-1242, 2004. 

J. Martinez and M. Silva. A simple and fast algorithm to obtain all invariants 
of a generalized Petri net. In C. Girault and W. Reisig, editors, Selected 
Papers from the First and the Second European Workshop on Application 
and Theory of Petri Nets, pages 301-310. Springer- Verlag, London, UK, 
1982. 

J. M. Mendez and R. Femat. Dynamic equivalence in tangent spaces from 
vector fields of chemical reaction networks. Chemical Engineering Science, 
pages -, 2011. submitted. 

A. L. Nagy, B. Szabo, T. Turanyi, and J. Toth. 201 1. in preparation. 

T. Nagy and Turanyi T. Uncertainty of Arrhenius parameters. International 
Journal of Chemical Kinetics, 43(7):359-378, July 201 1. 

D. Papp and B. Vizvari. Effective solution of linear Diophantine equation sys- 
tems with an application in chemistry. The Journal of Mathematical Chem- 
istry, 39(1):15-31, 2006. 

M. Rathinam, L. R. Petzold, Y. Cao, and D. T. Gillespie. Stiffness in stochastic 
chemically reacting systems: The implicit T-leaping method. The Journal of 
Physical Chemistry A, 119(24): 12784-12794, 2003. 

H. H. Robertson. The solution of a set of reaction rate equations, pages 178- 
182. Thompson Book Co., 1966. 

T. Sipos, J. Toth, and P. Erdi. Stochastic simulation of complex chemical re- 
actions by digital computer, I. the model. React. Kinet. Catal. Lett, 1(1): 
113-117, 1974a. 

T. Sipos, J. Toth, and P. Erdi. Stochastic simulation of complex chemical reac- 
tions by digital computer, II. applications. React. Kinet. Catal. Lett., 1(2): 
209-213, 1974b. 

G. Szederkenyi. Computing sparse and dense realizations of reaction kinetic 
systems. The Journal of Mathematical Chemistry, 47:551-568, 2010. 

G. Szederkenyi and K. M. Hangos. Finding complex balanced and detailed 
balanced realizations of chemical reaction networks. The Journal of Mathe- 
matical Chemistry, 49:1163-1179, 2011. 

J. Toth and V. Hars. Orthogonal transforms of the Lorenz- and Rbssler- 
equations. Physica 19D, pages 135-144, 1986a. 

J. Toth and V. Hars. Specification of oscillating chemical models starting form 
a given linearized form. Theor. Chim. Acta, 70:143-150, 1986b. 

J. Toth, G. Li, H. Rabitz, and A. S. Tomlin. The effect of lumping and expand- 
ing on kinetic differential equations. SIAM Journal of Applied Mathematics, 
57:1531-1556, 1997. 

J. Toth, A. L. Nagy, and D. Papp. ReactionKinetics — A Mathematica pack- 
age with applications i. Requirements for a reaction kinetics package. Chem- 
ical Engineering Science, 2011. this issue. 

J. Vardai and J. Toth. Hopf Bifurcation in the Brusselator. 

http://demonstrations.wolfram.com/HopfBifurcationInTheBrusselator/. 
from The Wolfram Demonstrations Project. 

A. I. Volpert. Differential equations on graphs. Mat. Sh., 88(130):578-588, 
1972. 



T. Wilhelm. The smallest chemical reaction system with bistability. BMC 
Systems Biology, 3(90):9 pages, 2009. 



